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I.  INTRODUCTION 


Computer  programs  for  calculating  the  half-space  admittance  and 
scattering  matrices  for  a finite  planar  rectangular  waveguide  antenna 
array  are  described  and  listed  in  this  report.  The  aperture  dimensions 
can  be  less  than  or  equal  to  the  feeding  waveguide  dimensions  and  the 
elements  are  uniformly  spaced  in  each  of  two  directions  in  either  a 
rectangular  or  isosceles  triangular  lattice.  The  general  theory  and 
method  of  computation  are  given  in  the  report  [1].  Equations  drawn 
from  reference  [1]  are  preceded  by  1-.  For  instance,  (1-5)  denotes 
equation  five  of  reference  [1]. 

In  summary,  the  procedure  is  an  application  of  the  method  of 
moments  to  an  integral  equation  formulation  of  the  problem.  The  un- 
knowns to  be  determined  are  the  coefficients  of  the  equivalent  magnetic 
currents,  which  are  proportional  to  the  tangential  electric  field  in  the 
aperture  regions.  A single  expansion  function  is  used  to  approximate  the 
electric  field  in  each  aperture. 

The  computer  program  subroutine  AY  which  calculates  the  modal 
TE  TM 

coefficients  A , A.,  (1-30,31)  and  the  characteristic  admittances 

TE  TM 

, Y^  (1-34,35)  is  described  and  listed  in  Section  II.  The  subroutine 

YHSP  which  calculates  elements  of  the  first  column  of  the  half-space 
tis 

admittance  matrix  Y (1-12)  is  described  and  listed  in  Section  III.  The 

subroutines  FCT  and  FCTA  which  supply  the  integrands  respectively  of  the 

single  integrals  of  (1-68)  and  the  double  integrals  of  (1-59)  to  the 

numerical  integration  subroutines  QG8  and  QG6A-QG6B  for  calculation  of 
hs 

Y are  described  and  listed  in  Section  TV.  The  subroutine  INT  which 

transfers  data  (limits  of  integrations  and  integrand  specification)  from 

subroutine  YHSP  to  subroutine  QG8  is  described  and  listed  in  Section  V. 

The  subroutines  QG6A-QG6B  and  QG8  which  are  respectively  a six-point 

double  and  an  eight-point  single  Gaussian  quadrature  numerical  integration 

subroutine  are  described  and  listed  in  Section  VI.  The  subroutine  CSMTZ 

wp  hs 

which  solves  equation  (1-15)  where  [Y  6 + Y 1 is  symmetric  Toeplitz 
(occurs  for  the  uniformly  spaced  linear  array  lattice  case)  is  described 
and  listed  in  Section  VII.  The  subroutines  MATMLT , TRMMLT,  MULTTR, 
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MATVCA,  and  LINSLV  which  solve  equation  (1-15)  where  [YWg  + Yhs]  is 
symmetric  block-Toeplitz  (occurs  for  the  uniformly  spaced  rectangular 
array  lattice  case)  is  described  and  listed  in  Section  VIII.  The 
subroutines  LINEQ  and  DECOMP-SOLVE  which  are  respectively  a complex 
matrix  inversion  and  a Gaussian  elimination  - LU  decomposition  complex 
linear  equation  solver  subroutine  are  described  and  listed  in  Section  IX. 
The  subroutine  TOPGEN  generates  a complete  matrix  (half-space  admittance 
or  scattering)  given  one  column  is  described  and  listed  in  Section  X.  A 
main  program  which  uses  subroutines  AY,  YHSP,  TOPGEN,  CSMTZ,  DECOMP, 
SOLVE,  and  LINSLV  to  obtain  the  half-space  admittance  and  scattering 
matrices  is  described  and  listed  in  Section  XI  along  with  sample  input- 
output  data. 


II.  DESCRIPTION  OF  THE  SUBROUTINE  AY 


The  subroutine  AY(A,Y)  stores  the  submatrices  defined  by  (1-30) 
and  (1-31)  in  A and  the  admittances  defined  by  (1-34)  and  (1-35)  in  Y, 

TE 

A^  is  stored  in  A(k  = (m+l)/2  + n/2  * LM) 

TF 

Y^  is  stored  in  Y(k  = (m+l)/2  + n/2  * LM) 

m = 1,3,5,...,  LM 
n = 0,2,4, ... , LN 
TM 

A^  is  stored  in  A(k  = LM*LN  + (m+l)/2  + (n-2)/2  * LM) 

TM 

Yk  is  stored  in  Y(k  = LM*LN  + (m+l)/2  + (n-2)/2  * LM) 

m = 1,3,5,...,  LM 
n = 2,4,6,...,  LN  . 


Note  that  only  odd  m and  even  n modes  are  calculated  and  that  the  first 
TE  TM 

subscript  of  A^k  and  A^k  in  (1-30)  and  (1-31)  has  been  dropped  since  the 

same  expansion  function  is  used  in  every  waveguide  region,  A = 

^2k  = •••  = = T^e  variables  AL,  AL1,  BL,  BL1,  LM,  LN,  ER,  and 

ETA  in  the  COMMON  statement  are  respectively  a/X,  a' A,  b/X,  b'/X,  LM,  LN, 

£ , and  r 
r 

given  by 


e , and  n where  X is  the  free  space  wavelength.  Minimum  allocations  are 


!! 
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L 


COMPLEX  Y(2  * LM  * LN  + LM) 

DIMENSION  A(2  * LM  * LN)  , F(LM) , STE(LN) , 
STM(LN) . 


DO  loop  10  stores 
b' 


2e 


a a ttA 


mrb ' 

sin  ■ 

, ,3  Vaa'  bb'  cos  2 nub’  in  STE  (2  X) 


2b 


^ V 

L'bTTX  V 


, nub' 

-z sin  — — - 

^ HIT  /b  • nrrrw 

aa ' bb^  COS  T "Sfb* ln  S™  (2} 


2b 


If  n * 0,  statement  30  calculates  STE(l)  using  sin(0)/0  = 1. 
stores 


DO  loop  12 


2a 


(Hi-  J-) 
v 2 .2’ 
a a 


, mTT  mTTa*  . _ ,111+1 N 

sin  T cos  -2J-  m F(— ) . 


If  |m  - a/a'|  < e (a  small  number),  statement  13  calculates  F((m+l)/2) 

2 2 2 2 
by  replacing  (l/(m  /a  - 1/a'  ))  cos(mua'/2a)  by  its  limit  (-  ua  / 4) 


TE  TM 

as  m approaches  a/a'.  DO  loop  20  calculates  the  coefficients  A,  , A.  , 
TE  TM  TE  k.  k 

, and  Y^  . Statement  31  stores  A^  in  A(k)  where 


A^E  = A(k)  = 


m 


+ <UT>2 


STE(|  + 1)  F(^) 


where 


k = (m+l)/2  + n/2  * LM  { 


f m = 1,3,5,...,  LM 
n = 0,2,4,...,  LN  . 


TM 


Statement  33  stores  A^  in  A(k)  where 


where 


= A(k) 


STE(^)  F(---— ) 

/m2  , n . 2 

V<2l>  + Zb 


k = LM  * LN  + (m+l)/2  + (n-2)/2  * LM 


1.3.5. . .. ,LM 

2.4.6. . ..,LN  . 


Statement  22  or  24  stores  Y of  (1-34)  in  Y(i)  where 


f" 

i = (m+l)/2  + n/2  * LM  -l 


1,3,5, ... ,LM 
0,2,4, .. . ,LN 


2 

If  the  calculated  value  of  (k./k)  - e is  zero,  then  statement  23 

2 -61  r TM 

replaces  (k^k)  - er  by  10  . Statement  32  stores  Y*  of  (1-35) 

in  Y(i)  where 


i = LM  * LN  + (m+1) /2  + (n-2)/2  * LM 


1.3.5. . .. ,LM 

2.4.6. . .. ,LN 


C LISTING  OF  THE  SUBROUTINE  AY 

C 

SUBROUTINE  AYCA.Yl 
COMPLEX  U.YI55J 

DIMENSION  A(50).F(5) .STE( 5) . STM I S) 

COMMON  AL1 • BL 1 *PI . PI  2, PI  3. U/ R 1 /ETA/R5/ AL . BL 
COMMON  /R6/ER.LM.LN.NT 
A I =SQRT  < 1 • / ( AL* AL 1*BL* BL 1 1 > 

A2=A  1*BL1.'<  AL*AL1  *PI  ) 

A3«1.4142t**A2 
A 4=1 ./<  AL*AL> 

AS* I ./( AL1*ALII 
A6*PI3*AL1/AL 
A 7=2.4  AL 
B1*BLI/BL 
82=B1*PI3 

B3=- 1.4 142 144 A 1*31/1  PI *AL1> 

84*2 .*BL 
ETA2=ETA*ETA 


I 

I 


■ 


1 


LN1=LN*1 

30  STE  ( 1 ) —A2 
JN1  = 2 

JN  2=2 

DO  10  JN=1.LN 
C 1=JN1*B2 

C2=S  IN ( C 1 )/Cl*<— l l**JM 
IF( JN.EQ.LNi  GO  TO  1 1 
STE<  JN2)=A3*C2 

11  STM* JN)=JN1 *B3*C2 
JN 1*  JN 1 ♦2 

JN2= JN2+ 1 
10  CONTINUE 
JM1  = I 

00  12  JM=1 • LM 
JM2= JM41 

IF1ABS* JM-AL/AH)  .LT.O.OOl)  GO  TO  13 
F(JM)=2./(JM1*JN1 *A4-A5)*C0S< JM1*A6>*< -1 >**JM2 
GO  TO  14 

13  FI JM )=— PI*AL*AL/2.*<-l l**JM2 

14  JMl=JMl+2 

12  CONTINUE 
KTE=0 

KTM=LN 1 *LM 
KTM1=LN*LM 
JN2=0 

OO  20  JN=1 • LN1 
JN1= JN-1 
BN= JN2/B4 
JM  1-  1 

OO  2 1 JM=1 . LN 
AM=JM1/A7 
C l=AM*AM >BN*8N 
C2=S0RT1C1 > 

C1=C  1-ER 
KTE=KTE4l 
IF(C1)  22.23,24 

22  Y(KTE ) =SQR T <— C 1 1 /ET  A 
GO  TO  25 

23  Cl=l.E-6*ER 

24  YCKTE)=-U*SQRT*C1 l/ETA 

25  1F(  JN.EQ.LNli  GO  TO  26 

31  A*KTE)=JM1*STE< JNJ*F< JM)/C2 

26  IF(JN.EQ.l)  GO  TO  27 
KT*=KTM* 1 

KTNl =KTM1 ♦ 1 

32  Y*KTMJ=EH/’<  Y*KTEi*ETA2  ) 

33  A1KTM1 )=STM< JNll*F( JM1/C2 

27  JM1=  JM1  42 
21  CONTINUE 

JN2=  JN2  4-2 
20  CONTINUE 
RETURN 
ENO 


_ _ _ , 

! 
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III.  DESCRIPTION  OF  THE  SUBROUTINE  YHSP 

The  subroutine  YHSP  (N,X,Y,YSP)  uses  the  single  integrals  of 

(1-68)  and  the  double  integrals  of  (1-59)  to  calculate  elements  of  the 

hs 

first  column  of  the  Y matrix  (1-12).  There  are  N apertures.  X(I) 
and  Y(I)  are  respectively  the  x and  y center  coordinates  of  the  Ith 
aperture. 

Minimum  allocations  are  given  by 

COMPLEX  YSP(N) 

DIMENSION  X(N) , Y (N)  . 

The  constants  D(l)  - D(6)  are  the  multiplying  factors  and 

K - K (1-61  to  64,  1-66  to  67)  of  the  integrands  of  (1-59), 

/ O 

D(l)  = X { ( 1 - l/4a'2)  (x.  - x.  + a’)  cos  -y  (x . - x.) 

l j a j l 

+ — (1  + l/4a'2)  sin  -y  (x . - x.)} 

tt  a'  j l 

D(2)  = X { ( 1 - l/4a'2)  (x.  - x.  + a')  sin  -y  (x . - x ) 

i 1 a J i 

- — (1  + l/4a'2)  cos  -np  (x.  - x.)} 

TT  a j 1 

D(3)  = (1  - l/4a'2)  cos  -y  (x.  - x ) 

a j i 

D(4)  = (1  - l/4a'2)  sin  -y  (x  - x.) 

a j l 

D(5)  = X { ( 1 - l/4a'2)  (x . - x + a')  cos  -y  (x.  - x.) 

1 i a j l 

- — (1  + l/4a'2)  sin  -y  (x . - x.)} 

tt  a j i 

D(6)  = X{ (1  - l/4a'2)  (x.  - x + a')  sin  -y  (x . - x ) 

J i a J 1 

+ — (1  + l/4a'2)  cos  -y  (x.  - x )} 
tt  a j i 


(a' , x . , 


and  x.  are  distances  per  unit  wavelength). 
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The  constants  TH1  - TH9  are  ngles  at  which  a rotating  line 
from  the  origin  (line  0-A  in  Fig.  1)  will  intersect  different  straight 
line  segment  junctions  defined  by  x = x_.  - x^  - a',  - x^,  x.  - x.  + a', 

and  y = Yj  - y±  - b\  y - yi>  y - y±  + b'  (see  Fig.  1). 

The  CALL  INT(N4,  N5,  N6,  PI,  P2,  XL,  XU)  statements  between  12 

and  11  select  the  correct  p (0)  and  p (0)  expressions  used  in  M,  (0)  - M,  (0) 

12  1 o 

(1-77  to  82),  select  the  integrand  for  each  subarea  (I-IV)  in  which  the 
single  integration  is  being  performed,  and  supply  the  0 limits  of  inte- 
gration. 


For 

N4  = 0, 
N4  = 1, 
N4  = 2, 
N5  = 0, 
N5  = 1, 
N5  = 2, 
N6  < 2, 

N6  = 2, 

N6  = 3, 


P1(6)  = 0 
p^(0)  = Pl/cos  0 
p^(0)  = Pl/sin  0 
P2(0)  = 0 
P2(0)  = P2/cos  0 
P2(0)  = P2/sin  0 


Integration  over  a subarea  I-IV  is 
still  being  performed 
Multiplying  factors  for  subareas  I 
and  III  are  used 

Multiplying  factors  for  subareas  II 
and  IV  are  used 


where  PI  and  P2  are  lines  which  define  the  subareas  of  integration  (I-IV) 

and,  therefore,  represent  the  limits  of  integration  for  the  variable  p. 

For  instance,  in  Fig.  1,  the  lines  which  define  the  limits  of  integration 

along  line  OA  in  subarea  IV  are  PI  = y.  - y.  and  P2  = y.  - y.  + b'.  XL  and 

j 1 J i 

XU  represent  respectively  the  0,  and  0 limits  of  integraf ~n. 

lower  upper 

The  arguments  of  the  CALL  INT  statements  will  vary  since  the  area 
of  integration  (I-IV)  can  move  throughout  the  first  quadrant  of  the  uv 
plane  (assumes  x^  _<  x^).  There  will  be  overlapping  into  other  quadrants 
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when  = x^,  y ± ? Yy  xt  * x..  , y.^  = y^  , and  x±  = , y±  = y^  . There 

are  four  specif ic  cases  to  consider  for  determining  the  integrations  to 
be  performed: 

Case  1)  x^  f x . , y^  f y^  - this  integration  will  require  twelve 
CALL  INT  statements  representing  the  twelve  integration 
areas  of  subareas  (I-IV) . The  twelve  statements  are 
located  between  statements  13  and  14. 


Case  2) 


Xi  = Xj’  yi  ^ ~ this  integration  will  require  only 

four  CALL  INT  statements  because  of  symmetry  con- 
siderations. The  four  statements  are  located  between 
statements  14  and  15. 


Case  3) 


Xi  ^ X j ’ yi  = yj  ~ integrati°n  like  case  2 will 

require  only  four  statements  which  are  located  between 
statements  15  and  11. 


Case  4) 


Xi  = Xj ’ yi  = yj  ~ this  integration  which  represents 
the  calculation  of  the  self-admittance  will  require 
only  two  CALL  INT  statements  which  are  located  between 
statements  12  and  13. 


If  (Xj  - x^)  > 4a,  a solution  involving  the  double  integrals 
of  (1-59)  is  used.  The  calling  statement  is  CALL  QG6A(XL,XU,YL,YU,N9,Z) . 
The  variables  XL,  XU,  YL,  YU  represent  respectively  the  lower  and  upper 
limits  of  integration  for  the  u and  v variables  of  integration  in  (1-59) . 


If  N9  = 1,  the  first  of  the  four  double  integrals  of  (1-59)  are 
evaluated  (if  N9  = 2,  the  second,  etc.).  There  are  three  specific  cases 
to  consider  for  determining  the  integrations  to  be  performed: 

Case  1)  x^  5s  x . , y ^ y . - this  integration  will  require 
four  CALL  QG6A  statements  representing  the  four 
subareas  of  integration  (I-IV) . The  four  state- 
ments are  located  between  statements  11  and  25. 
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Case  2) 


x.  = x.,  yi  ^ y.  - this  integration  will  require  only 

two  CALL  QG6A  statements  because  of  symmetry  considerations. 

The  two  statements  are  located  between  statements  25  and  26. 


Case  3)  ^ x^,  y^  = y^  - this  integration  like  case  2 will 

require  only  two  CALL  QG6A  statements  which  are  located 
between  statements  26  and  10. 


C LISTING  OF  THE  SUBROUTINE  YHSP 

C 

SUBROUTINE  YHSP ( N .X. V. YSP » 

COMPLEX  G3.G5.G7.G9.SI 1 0 ) • SA • SB • U. UI .YSPI6) 
DIMENSION  X(6)«Y<6> 

COMMCN  A1.B1 »PI«PI2»PI 3»  U/Rl /ET  A/R2/DI 6 I • S 
COMMON  /R3/G 1 • G2 • G3  » G4  «G 5«  G6 »G7. G8*G9*  G i 0 • G1 1 
COMMON  /R5/AL.BL/R7/YY3. YY4 
Ul=U/< A1*B1*ETA) 

OO  10  J= 1 • N 
X1  = X(I  ) 

X2=X< J> 

Y1=Y  (1  ) 

Y2=Y( J> 

XX1=X2-XI-A1 
XX2=X2-X1 
XX3=X2-X1  + A1 
XX4=X1-X2*-A  1 
YY 1= V2— Y 1 — B 1 
VY2=  Y2— Y 1 
YY3=Y2-YI+B  1 
YY4=Y1-Y2+B1 
XA=(  X2-XD/G1 
DI=4.*AL 
Cl=COS« XA» 

C2=S INI XA) 

D«3)=C1*G10 

D14)=C2*G10 

C1=C1*G1 1 

C2=C2*Gl 1 

Of  1 )=XX4*0( 31+C2 

DC2)=XX4*DC4)-C1 

D I 5)=XX3*D( 31-C2 

016) =XX3*0 ( 4 ) +C1 

IFIXX2.GT.DI)  GO  TO  11 

TH1  = AT  AN2I YY1 • XX3 I 

TH2=  AT  AN2I Y Y2 • XX3 ) 

TH3= AT  AN2I YY3.XX3) 

TH4-AIAN2( YY1.XX2) 

TH6=ATAN2< YY3.XX2 ) 

TH7*ATAN2<  VY1.XX1 ) 


11 


12 


13 


16 


17 


18 


19 


20 


21 


22 


23 
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T HQ=  AT  AN2 1 Y Y2 . XX 1 ) 

TH9=ATAN2( Y Y3 • XXI 1 
IF1J.EQ.1)  GO  TO  12 
THS=ATAN2< YY2.XX2) 

IF1J.GT.1)  GO  TO  13 
CALL  INT ( 0 • 1.0.0. . XX3. O.  » TH3) 
CALL  INT<0 .2.3.0. .YY3.TH3.PI3) 
YSP< J)=4.*( YY3*St 9)-S( 10)1 
GO  TO  10 

IF!  1 Xl.EQ.X2). AND . 1 Y 1 . NE . Y2  ) ) GO 
IF! 1X1.NE.X2).AND.1Y1.EQ.Y2) ) GO 
IF  1 T H7 . GT  . T H5 1 GO  TO  16 
CALL  INTO. 1 ,0,YV 1.XX2 .TH4.TH7) 
CALL  INTO  . 1. 1. XXI. XX2.TH7.TH5) 
CALL  INTO  .2.2.XXI.YY2.TH5.TH8) 
GO  TO  17 

CALL  INTO. 1.0.YYI.XX2.TH4.TH5I 
CALL  INT12.2.1.YYI.YY2.TH5.TH7) 
CALL  INK  1 .2.2.XX1.YY2.TH7.TH8) 
SA<=YY4*S19) 

SB=S (10) 

IF1TH8.GT.TH6)  GO  TO  18 
CALL  INTI2# 1*0. YY  2 • X X2 .TH5.TH8) 
CALL  INTI 1 • 1 .1 .XX1.XX2.TH8. TH6) 
CALL  INTO  .2.2.XX1.YY3.TH6.TH9) 
GO  TO  19 

CALL  I NT (2.1 . 0 . YY 2. XX2 • T H5»  TH6 ) 
CALL  I NT (2.2.1. YY2. YY3 . TH6»  TH8  J 
CALL  INK  1 .2.2.XX1.YY3.TH8.  TH9) 
SA=SA*YY3*S19) 

SB=  S B—  S(  10) 

YSP( J ) =SA+SB 

IF ( TH4.GT. T H2)  GO  TO  20 

CALL  1 NT  < 2 • 1.0. YY I • XX3.TH I , TH4 ) 

CALL  I NTO  . 1 ,1  .XX2.XX3.TH4.TH2) 

CALL  INTI  1 .2.3. XX2.YY2.TH2.THS) 

GO  TO  21 

CALL  INTI2.1.0.YY1.XX3.TH1.TH2) 
CALL  I NT (2.2* l.YY 1.YY2. TH2.  TH4) 
CALL  INTO  .2*3.  XX2.  YY2  . T H4.  TH5  ) 
SA=Y  Y4*  S( 9 ) 

SB=SI10) 

IFITH5.GT.TH3)  GO  TO  22 
CALL  INTI2.1.0*  YY  2. XX3 .TH2.TH5) 
CALL  I NT (1.1.1 . XX2 . XX3 .TH5.TH3) 
CALL  INK  1.2.3. XX2.YY3.TH3.TH6) 
GO  TO  23 

CALL  (NTI2.1.0.YY2.XX3.TH2.TH3) 
CALL  I NT  12.2.1. YY2.YY3. TH3. TH5 ) 
CALL  I NT (1.2.3. XX 2. YY3 • T H5»  TH6 ) 
SA=SA4-YY3*S19) 

SB=SB— SI  10) 

YSPl J)=YSP1 J)+SA»SB 
GO  TO  10 

CALL  INK2.2.0.YY1.YY2.PI3.TH8) 
CALL  INT12.1 .2.YY1.XX1 .TH8.TH7) 
SA»YV4*SI 9) 

SB=S (10) 


TO 

TO 


14 

15 


i 
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CALL  INTC2.2.0.YY2.YY3.PI3.TH9) 

CALL  INTC2, 1 .2.YY2.XX1 ,TH9. TH8J 
SA=SA4YY3*SC9» 

SB=SB-S( 10) 

YSPC  J)=2.*C  SA  + SBi 
GO  TO  10 

15  CALL  INT(1 .2.0. XXI. YY1 .TH7.TH4) 

CALL  I NT  11*1 *2,XX I • XX2 • T H4»  0 • ) 

VSPI  J)=yV4*S19)»SU0) 

CALL  INTC 1.2,0.XX2.YYt .TN4.TH1) 

CALL  INTO  .1  .3.XX2.XX3.TH1.0.) 

YSPC  J)  = 2.*C  YSPC  Ji ♦YY4*S<  9)* SC  10)  ) 

GO  TO  10 

II  IF<< Xl.EQ.X2) .ANO.C Y1.NE.Y2) ) GO  TO  25 
IF< C XI .NE.X2) .ANO.C Yl. EQ.Y2) ) GO  TO  26 
CALL  QG6A (XXI  * X X2  . Y Y 1 * YY  2 • 1 « SA I 
SB=SA 

CALL  QG6  AC  X X2. XX3  » YY 1.YY2.2* SA ) 
SB=SB*SA 

CALL  QG6A(XX1,XX2.YY2.YY3.3.SA> 
SB=SB+SA 

CALL  QG6AC  X X2» XX3 .YY2.YY3.4*  SA) 

YSPC  J>=SAC-SB 
GO  TO  10 

25  CALL  0G6AC XXI .XX2.YY1, YY2. 1 . SAI 
SB=2 •*  SA 

CALL  QG6AC XXI. XX2 . YY2 . YY3 , 3 . SA » 

YSPC  J)=2.*SA«-SB 
GO  TO  10 

26  CALL  QG6AC  XXI  * XX2 • YY1»  YY2, 1 • SA ) 
SB=2.*SA 

CALL  QG6AC XX2. XX3.YY1,  YY2.2 .SAI 
YSPC J)=2.*SA*SB 
10  CONTINUE 

OO  30  1=1, N 
YSPC  1 >=U1*YSPC 1 } 

30  CONTINUE 
RETURN 
END 
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IV.  DESCRIPTION  OF  THE  SUBROUTINES  FCT  AND  FCTA 

The  subroutine  FCT  (W,X)  supplies  the  integrands  for  the  single 
numerical  integration  subroutine  QG8 . The  W(I)  are  the  integrands  of 
equation  (1-68)  evaluated  at  0 = X where  X is  the  variable  of  integra- 
tion for  QG8. 

The  logic  between  statements  40  and  6 determines  the  p^(0)  and 
P2(0)  expressions  used  in  M^©)  - M^(d)  ( 1-77  to  82). 

For 


N4 

= o , 

C = 

p1(6) 

= 0 

N4 

= 1 , 

C = 

P1(0) 

= Pl/cos 

e 

N4 

= 2 , 

C = 

P1(0) 

= Pl/sin 

0 

N5 

= o , 

D = 

P2(0) 

= 0 

N5 

= 1 , 

D = 

P2(0) 

= P2/cos 

0 

N'5 

= 2, 

D = 

P,(0) 

= P2/sin 

e 

W(l)  - W(8)  represent  respectively  L..  (0)  - Lo(0)  (1-69  to  76). 

1 O 

The  logic  statements  between  41  and  37  supply  the  limiting  expres- 
sions (1-84  to  99)  when  l-^p  cos  0 + k|  < e (a  small  value).  For  a 
v ^ 

first  order  zero  (—  cos  9 + k)  in  the  denominator  of  L. (0)  - Lo(0) 

3 1 o 

(1-69  to  76),  the  e value  of  El  = 0.0002  was  used  while  for  higher  order 
zeros,  e = E2  = 0.002  was  used. 

The  subroutine  FCTA(N9,X,Y,Z)  supplies  the  integrands  for  the 
double  numerical  integration  subroutines  QG6A  and  QG6B . N9  specifies 
which  of  the  four  double  integrals  of  (1-39)  will  be  evaluated  (if 
N9  = 1,  the  first  of  the  four,  if  N9  = 2,  the  second,  etc.).  X and  Y 
are  the  variables  of  integration  respectively  for  QG6A  and  QG6B.  Z 
is  the  integrand  of  one  of  the  four  double  integrals  of  (1-59)  evaluated 
at  X and  Y.  The  integrands  of  the  double  integrals  of  (1-59)  are 
transferred  from  subroutine  YHSP  to  subroutine  FCTA  through  the  C0MM0N/R2/ 
and  C0MM0N/R7/  statements. 


u o 


LISTINGS  OF  THE  SUBROUTINES  FCT  AND  FCT  A 

14 

SUBROUTINE  FCTCW. X» 

COMPLEX  E(8).F(8i .G3 .G5 • G7 • G9 .SA • SB . U. W <8 i 

COMMON  A1.B1 .PI.PI2.PI3*U/R3/G1.62.G3.G4.G5.G6.G7,G8,G9*G1 0«GII 

COMMON  /R4/N4 .N5.N6.Pl .P2.XL.XU 

Cl=COSCX> 

C2— S IN ( X ) 

E 1 =0  *0002 
E2=0  *002 
P=CI /Al 
S=PI*CP-2. ) 

T=Pl *CP+2. ) 

40  IF I N4. NE . 0 ) GO  TO  1 
C=0. 

GO  TO  3 

I  IFIN4.E0.2)  GO  TO  2 
C=P1/C1 


GO  TO  3 

2 C=PI/C2 

3 I F ( N5*  NE . 0 ) GO  TO  4 

0=0. 

GO  TO  6 

4 IF ( NS. EQ.2 ) GO  TO  S 
D=P2/CI 

GO  TO  6 

5 0=P2/C2 

6 Y=CI-2.*A1 
P=D*S 

EC l J=COSCP>  +U*SIN«P) 

P=C*S 

EC2)=COSCP> +U*S  INC  P) 

Z=Cl  +2.*A1 
P=D*T 

FC 1 > =COSC  P >- U* SI N CP) 

P=C*T 

FC2)=COSCPI-U*SINCP» 

EC  3) =C  EC  I ) — EC  2 ) >/V 
FC3)=CFC1 >-F<2> )/Z 
SA=D*E  C 1 ) 

SB=C  *E  C 2 ) 

EC4)=CSA-SB)ZY 

SA=0*SA 

SB=C4SB 

EC6)=CSA-SB»/Y 
SA=0*FC I ) 

SB=C  *F C 2 I 
FC4»=C  SA-SBI /Z 
EC  5)=EI 3 )/Y 
FC5»=FC3)/Z 
SA=0*S A 
SB=C*SB 

FC6»=C  SA-SBI/Z 
EC  7I=EC  41/Y 
FC7)=FC4)/Z 
EC8) =EC  5JZY 
FC8»=FC  5)/Z 

41  IFCABSCYI.LT.Ell  GO  TO  20 
IFCABSCZ1.LT.E1)  GO  TO  21 
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GO  TO  22 

20  WC1)=-G3*FC3»+0.54CO-C> 

GO  TO  22 

21  WC  1 MG3*EC3)«-0.5*CD-C» 

22  IFCABSCYI.LT.E1)  GO  TO  23 
IFCABSCZI.LT.El)  GO  TO  24 
W (2  1=—  G2*C  El  3 1 4FC  3)  ) 

GO  TO  25 

23  MC 2)=— G2*F ( 31— 0«5*U* CO— C ) 

GO  TO  25 

24  WC  2 ) =— G2*E 13) *0«5*U*  CD— Cl 

25  IFC ABSC Y1.LT.E2)  GO  TO  26 
IFCABSC Z1.LT.E2)  GO  TO  27 
tfCSl=G3*CEC  41-FC4)  )+G6*CEC5)*FC5  )l 
GO  TO  26 

26  H=0.25*CD*D-C«C1 

W( 5 } =— G3*F  C 4 1 ♦G6*F( 5>*H 
GO  TO  26 

27  H=0. 25*CD*D-C*C1 

M ( 5)=G3*E( 4 1 *G6*E (5) FH 

28  WC31=CI*WC5) 

MCS1=C2*»C 5) 

IFCABSCY1.LT.E21  GO  TO  29 
IFCABSCZ).LT.E2)  GO  TO  30 
W ( 61*=—  G2*(  E C 4 1 4F  C 4)  M-G7*CEC  Sl-FC  51  > 

GO  TO  31 

29  MC6)=-G24FC 4 1-G74FC 5 1- U*H 
GO  TO  31 

30  WC6)=-G2*EC41*G7*EC5)*U»H 

31  WC4)=C1*WC6) 

VC61=C2*W<61 
C1=C  1*C2 

IF(ABS( Yl.LT.E21  GO  TO  32 
IFCABSCZ1.LT.E2)  GO  TO  33 

If  C 71 =C 1*<G3*CE( 6 1— FC 6)  1»G4*CEC71 ♦FI  7) ) -G9* CE 1 8 » -F < 8 ) 1 ) 

GO  TO  34 

32  H=CD*D*D-C*C*Cl/6. 

**C71=Cl*C-G3*FC6i*G4*FC  71*G9*FC  814H1 
GO  TO  34 

33  H=(D*D*D-C*C*Cl/6. 

«C7)=Cl*CG3*EC61+G4*E( 7}-G9*E< 81 +H1 

34  IFCABSC Y1.LT.E2I  GO  TO  35 
IFCABSCZ1.LT.E21  GO  TO  36 

WC8  >=C1 *C-G2*CEC61*FC6 I > 4G5*C  EC  7 )-F<  7 ) ) ♦G8*C  EC  8) *FC  81)  1 
GO  TO  37 

35  W C 8)=CI *C-G2*FC61-G5*FC  71 *G8*FC 8 1-U4H1 
GO  TO  37 

36  «C8)=C1*C-G2*EC6>»G5*EC  71  *G8*E  C 8 1 «-U*H  1 

37  CONT  INUE 
RETURN 
ENO 

SUBROUTINE  FCT AC N9. X . Y . Z 1 
COMPLEX  G.G3.G5.G7.G9. SC  10) .U.Z 
COMMON  AL1 . BL 1 .P I .P 12. P I 3. U/R2/D C6 1 • S 

COMMCN  /R3/G1 • G2. G3. G4 • G5 . G6 . G7. G8.G9.  G10. G1 1/R7/YY3.YY4 
P 1=X/G1 
ClcCOSCPl ) 

C2*SINC  Pll 

C3=SQRT(X*X*Y*Vl 

C4=PI2*C3 


G=tC0S<C4)-U*SIN< C4>  >/C3 
IFIN9.EG.2)  GO  TO  1 
IFIN9.EQ.3)  GO  TO  2 
IF ( N9 . EQ • 4 ) GO  TO  3 

Z=G*CYY4+Y) *<CI*<  0(  1 l+X*0(3J  J «-C2*<  D<  2 ) +X*0<4)1J 
RETURN 

1 Z=G*(YY4*Y)*(Cl*IO(5>-X*Dl3l ) +C2 *< D< 6 ) - X*D ( 4 ) » ) 
RETURN 

2 Z=G*<YV3-Y  ) *<C  1*1  Dt  l )+X*D(3)  I«-C2*<  D(2)«-X*D(4))I 
RETURN 

3 Z=G*< YY3-Yi*(C 1*1 01 5 )-X*0C3»  ) *C2* < OC 6 > -X4D ( 4 ) » ) 
RETURN 

END 


n o 
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V.  DESCRIPTION  OF  THE  SUBROUTINE  INT 

The  subroutine  INT(N4,  N5,  N6,  PI,  P2,  XL,  XU)  transfers  the  data 
N4,  N5,  N6  which  is  supplied  by  the  subroutine  YHSP  and  required  by  the 
subroutine  FCT  through  the  C0MM0N/R4/  statement,  calls  the  single 
numerical  integration  subroutine  QG8,  and  calculates  the  S ( 9)  and  S(10) 
terms  which  represent  the  four  single  integrals  of  (1-68) . The  arguments 
N4,  N5,  N6,  PI,  P2,  XL,  and  XU  are  defined  in  the  description  of  the  sub- 
routine YHSP  and  will  not  be  repeated  here. 

1 

The  following  logical  statements  are  located  between  statements 
3 and  4.  If  N6  is  less  than  2,  no  S(9)  or  S(10)  term  is  calculated  (one 
or  two  more  intermediate  integrations  are  required  before  evaluation) . 

Also,  if  N6  = 0,  the  variable  S which  represents  the  result  of  inter- 
mediate integrand  evaluations  of  subroutine  QG8  is  reset  to  zero.  If 
N6  = 2,  S ( 9 ) represents  the  first  single  integral  (I)  of  (1-68)  and  S(10) 
represents  the  third  (III).  If  N6  = 3,  S(9)  represents  the  second  inte- 
gral (IT)  of  (1-68)  and  S(10)  represents  the  fourth  (IV). 


LISTING  OF  THE  SUBROUTINE  INT 

SUBROUTINE  INTIN4.N5.N6.P1 .P2.XL.XU) 

COMPLEX  S(  10).T(8) 

COMMON  /R2/0I6)*S/R4/N7.N8,N9.P3«P4.XL1 .XU I 

N7=N4 

N8=N5 

N9=N6 

P3=P  1 

P4=P2 

XL!  = XL 

XUi  = XU 

IFIN6.NE.0)  GO  TO  I 
OO  2 1*1,8 
SI I > = ( O.  .0.  ) 

2 CONTINUE 

1 CALL  QG8I T ) 

OO  3 1*1 .6 
SI  1 1 *S< I I ♦ T I I ) 

3 CONTINUE 
IFIN6.LT.2)  GO  TO  4 

IFIN6.EQ.2)  SI  91*01  ll*SIl)*DI2)*St2)FDt3»4SI3>4DI4)*SI4» 
IFIN6.E0.3)  SI9)*DI5J*S<  1 )40 I 6 1 4SI 2 l-DI 3) * SI  3) -O I 4 ) * SI  4 ) 
IFIN6.E0.2)  SI 10)=O<  I » *S  I 5 > »D  1 2 > *S I 6 » *0 1 3 > *S  I 7 » *OI  4 > • SI  8 > 
IFIN6.EQ.3)  SI  101*01 51 *S( 5)*OI6l •SI6I-DI3 )*SI 7 I — 014) *SI8 ) 


4 CONT INUE 
RETURN 
END 
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VI.  DESCRIPTION  OF  THE  SUBROUTINES  QG6A,  QG6B,  AND  QG8 

The  combination  of  subroutines  QG6A(XL,  XU,  YL,  YU,  N9,  Z)  and 
QG6B(X,  YL,  YU,  N9,  Z)  perform  a six-point  Gaussian  quadrature  double 
numerical  integration  on  the  double  integrals  of  (1-59) . The  basic 
integral  under  consideration  can  be  represented  as 


XU 


YU 


I = 


dx 


dy  F(x,y) 


(1) 


XL 


YL 


This  integral  can  be  transformed  from  the  integration  intervals 
(XL,  XU),  (YL,  YU)  to  the  intervals  (-1,  1),  (-1,  1)  by  the  trans- 
formation 


so  that  (1)  becomes 


i - (3fc”=)  (Sts, 


-1 


2x  - (XU  + XL) 

(2) 

A 

(XU  - XL) 

. 2v  - (YU  + YL) 

(3) 

y 

(YU  - YL) 

1 

, r XU-XL  XU+XL. 

rYU-YL  , , YU+YL 

,t  2 y + 2 

dx' 

dy'  F([--2-  --  x'+  2 ] 

-] 

L 

(4) 

A double  six-point  Gaussian  quadrature  formula  [2,  Section  7.2]  is  used 
to  evaluate  (4) , 

I - (*!=£, (IML)  l l A<6>A<«  F(x<6>,  y<6)>  (5) 

z i=l  j=l  J J 

where  x^6\  are  the  roots  of  the  Legendre  polynomial  of  degree 

6 _ P (x<6>  or  y^)  = 0 (the  Legendre  polynomials  are  orthogonal  on 

6 1 j (6)  (6)  . CJ  . 

the  (-1,  1)  interval).  The  coefficients  A.  , A are  defined  as 


.(6) 


,(6) 


3 P'  (xf')  P5(x‘6)) 


3 P6  (yj6>)  P5(yj6)) 


(6) 


(7) 
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('  represents  derivative)  and  are  tabulated  [2,  p.  337]. 

In  the  calling  arguments,  (XL,  XU,  YL,  YU)  represents  respec- 
tively the  lower  and  upper  limits  of  integration  for  the  x and  y 
variables,  N9  is  used  in  the  subroutine  FCTA  to  specify  the  function 
to  be  integrated  (if  N9  = 1,  the  first  double  integral  integrand 
of  (1-59)  is  specified,  if  N9  = 2,  the  second,  etc.)  and  Z is  the 
result  of  the  integration. 

Subroutine  QG8(Y)  performs  a single  eight-point  Gaussian 
quadrature  numerical  integration  on  the  single  integrals  of  (1-68) . 
In  the  calling  argument,  Y is  the  result  of  the  integration.  The 
variables  N4,  N5,  N6,  PI,  P2,  XL,  and  XU  in  the  C0MM0N/R4/  statement 
are  defined  in  the  description  of  subroutine  YHSP  and  will  not  be 
repeated  here. 


C LISTINGS  OF  THE  SUBROUTINE  QG6A , QG6B , AND  QG8 

C 

SUBROUTINE  QG6A< XL . X U. YL . YU . N9 . Z > 

COMPLEX  Z.Z1.Z2 
A=0.5*<  XL-*- XU) 

B=XU— XL 

C=0. 4t>62348*B 

X=A+C 

CALL  QGEBt X.YL.YU.N9.Z1 ) 

X=  A — C 

CALL  QGbBtX .YL.YU.N9.Z2) 

Z=0. 08S66225* ( Zl* Z2> 

C=0. 330604 7 *B 
X=A+C 

CALL  QG6BI  X , YL  , YU  #N9  ,Z  1 ) 

X=  A — C 

CALL  QG6BI X ,YL  ,YU.N9 .Z2) 

Z=Z+0. 1 6038 08*1 Zl +Z2) 

C=0. 1 1 93096  *B 
X=A-*C 

CALL  QG6bCX,YL.YU.N9.Zl) 

X=A-C 

CALL  QG6B< X.YL.YU.N9.Z2) 

Z=B* < Z*0 #2339570* ( Z I *Z2) ) 

RETURN 

END 

SUBROUTINE  QG6B(  X»YL»YU»  N9»  Z ) 

COMPLEX  Z • Z 1 • Z2 

A=0. 5*1 VL+YUJ 

B=YU-YL 

C=0 • 466234  8 *B 

Y = A4-C 

CALL  FCTA<N9.X,Y.Z1 ) 
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Y=A-C 

CALL  FCTA(N9,X.Y,Z2) 

Z=0. 0856622 5* ( Z1FZ2) 

C=0 • 330 604  7 *B 
Y=A  + C 

CALL  FCTA(N9.X.Y.Z1) 

Y-A-C 

CALL  FCTA(N9.X.Y,Z2> 

Z=ZF0. 1803808*1 Z1FZ2) 

C=0. 1 X 93096*0 
Y=AFC 

CALL  FCTA(N9.X.Y.  Z1  ) 

Y=A-C 

CALL  FCTA< N9.X.Y. Z2) 
Z=B*(ZF0.2339570*(Z1 FZ2)  > 

RETURN 

END 

SUBROUTINE  QG8 ( Y ) 

COMPLEX  T(8) • T 1 ( 8 ) » Y ( 8 ) 

COMMON  /R4/N4.N5.N6.P1 .P2.XL.XU 
A=.5*< XL+XU) 

B=XU— XL 

C=.4801449*B 

X=A+C 

CALL  FCT(T.X) 

DO  1 1=1 .8 

Till  > = T (I  ) 

1 CONTINUE 
X=  A— C 

CALL  FCT(T.X) 

DO  2 1=1.8 

Y(I  )=.05061427*<T< I) FT  1(1)) 

2 CONTINUE 

C=. 3983332 *B 
X=AFC 

CALL  FCT(T,X) 

DO  3 1=1.8 
TUI  ) = T<  I ) 

3 CONT INUE 
X= A-C 

CALL  FCT(T.X) 

DO  4 1=1.8 

Y ( I ) =Y (I)F.lll 1 905* ( T(  I ) FT1 ( I ) ) 
♦ CONT INUE 

C=.2  62  7662*B 
X=A*C 

CALL  FCT(T.X) 

DO  5 1=1.8 
Till )=T ( I ) 

5 CONT INUE 
X=A-C 

CALL  FCT(T.X) 

DO  6 1=1 .8 

Y ( 1 ) =Y ( I ) F • 1 568533* ( T<  I ) FT  1(1)) 

6 CONTINUE 
C=.091 7 1 732*B 
X=  AFC 

CALL  FCT(T.X) 
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DO  7 Isl.e 
T 1 ( X )=T<  a 

7 CONTINUE 
X=  A-  C 

CALL  FCT(T.X) 

DO  8 1=1.8 

Y(I)=8*(Y(I)*.I8I3419*<T(I>*TI(I))) 

8 CONTINUE 
RETURN 
END 


' 
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VII.  DESCRIPTION  OF  THE  SUBROUTINE  CSMTZ 

The  subroutine  CSMTZ(N,  A,  B,  X)  solves  the  set  of  equations 


lnsn 


(8) 


where  is  an  N * N complex  symmetric  Toeplitz  matrix.  In  the  argument 

of  CSMTZ,  N is  the  number  of  unknowns  in  (8),  A is  the  input  array  whose 

elements  are  the  first  row  of  L^,  B is  the  input  array  whose  elements  are 

those  of  dxl,  and  X is  the  output  array  whose  elements  are  those  of  s.,.  It 
N N 

is  assumed  that  (8)  is  normalized  so  that  the  first  element  of  the  first 
row  of  is  equal  to  unity. 

This  subroutine  solves  equation  (1-15) 

[YWg  + Yhs]V  = Ilmp  (1-15) 

where  [YWg  + Y^s]  is  a symmetric  Toeplitz  admittance  matrix  (linear 
array  lattice  case)  . 

The  two  main  advantages  of  the  algorithm  used  in  this  subroutine 

are 


(1)  It  solves  equation  (1-15)  directly  without  inversion 
of  [YWg  + Y^S]  requiring  roughly  2N^  multiplications 
and  divisions. 

wc  hs 

(2)  It  requires  only  the  first  now  of  [Y  + Y ] and,  there- 
fore, minimizes  the  storage  requirement  for  the  subroutine. 

The  derivation  of  the  algorithm  can  be  found  in  [3,  4]  and  will 
not  be  repeated  here.  However,  the  algorithm  logic  will  be  presented 
using  the  notation  developed  in  [3,  4].  The  following  notation  will 
be  used  (the  same  used  by  Zohar  in  [3,  4]):  Greek  letters  are  used  for 
scalars,  capital  letters  for  square  matrices,  lower  case  letters  for 
column  matrices,  ~ denotes  transpose,  and  ^ is  a reversal  symbol  - g^ 
denotes  the  reversed  order  of  g^»  that  is,  = (gjc)jt+j_^  j • 

The  matrix  of  (8)  is  bordered  as  follows. 
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L 


S 


N-l 


u N— 1 


N-l 


(9) 


where  r„  , = [A(2) , A(3),...,  A(N) ] (A(i)  Indicates  the  ith  component 
N-l 


(numbered  columnwise)  of  array  A).  In  (8),  d^  = [6^,62> . . . ,6^] . The 


algorithm  is  based  on  a recursion  relation  with  initial  values  given  by 


S1  = 61  ’ pl  = ~ A(2)  ’ \ = 1 ~ A(2)  * A(2)  • (10) 


Recursion  of  s^,  e^,  and  X^  for  i-1 , 2 , . . . ,N-2  is  given  by 


9i  = 6i+l  " §i  fi 


(11) 


D.  = - A(i+2)  - ? e. 
1 11 


(12) 


Si+1 


Si  + (9,/X.)  e. 


8,/X. 

l i 


(13) 


'i+1 


ei  + (W  Si 


VAi 


(14) 


Xi+1  = Xi  - ni/Xi 


(15) 


The  last  computed  values  are  ^-2’  SN’  eN-l’  anc^  XN-1‘  Note 


that  a^,  g^,  and  y appearing  in  [3]  are  not  needed  because  of 


symmetry. 

Minimum  allocations  are  given  by 

COMPLEX  A(N) , B(N) , E(N) , ES (N-l) , X(N) . 
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The  initial  values  for  recursion  are  first  computed:  s^  = X(l),  = E(l)  , 

and  A^  = LA.  DO  loop  10  is  the  main  recursion  loop  whose  index  is 
i(i=l,2, . . . ,N-2) . DO  loop  11  computes  SI  = s^  and  El  = ? e^.  Following 

statement  11,  and  0^  are  computed.  DO  loop  12  computes  s.+^  and  stores 
e^  in  dummy  array  ES  because  e^  is  needed  in  DO  loop  13  to  recompute  itself. 
DO  loop  13  computes  e^.  Following  13,  A^+^  is  computed  and  the  index  i of 
DO  loop  10  is  stopped.  DO  loops  14  and  15  compute  s^. 


C LISTING  OF  THE  SUBROUTINE  CSMTZ 

C 

SUBROUTINE  CSMTZ I N, A *B * X > 

COMPLEX  A(6 ).B(6). EI6I.ES (5), X(6» 
COMPLEX  El .ET.ET1 .LA ,S 1 .TH.TH1 
XII )=B( 1 ) 

E(1 ) =— A ( 2 ) 

LA=1  •—  A < 2 ) *A (2) 

Nl=N-2 

N2=N-1 

DO  10  I = 1 » N 1 
El=( 0. .0.) 

S 1= I 0..0.) 

DO  11  J=  1 . I 

S1  = SH-X(  J1  *AI I-JF2J 

E1=E l+EI I-J+l )*A{ J+l > 

11  CONTINUE 
ET=-A< I+2I-E1 
TH=BCI  + 1 ) — S 1 
TH1=TH/LA 
ET1=ET/LA 

DO  12  K= 1 . I 

X(K)=X<K)*TH1*E< I-KFl) 

ES(K  ) = E ( K ) 

12  CONTINUE 
DO  13  K= 1 . I 

ECK)=ES(KM-ETI  *ES<  I-KF  1 ) 

13  CONTINUE 

X ( I + 1 ) = T H/L  A 
E(  14-1  )=ET1 
LA=LA-ET»ET 1 
10  CONT INUE 
Sl=< 0. .0. ) 

DO  14  J=1,N2 

S1=S  14-X  ( J)  * A<  N-J4-U 

14  CONTINUE 
TK=B(N1-S1 
TH1  = TH/LA 

DO  15  K=t.N2 
X(K)=XCK)+TH1*E(N-K) 

15  CONTINUE 
XIN2+1 ) = TH1 
RETURN 

END 
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VIII.  DESCRIPTION  OF  THE  SUBROUTINES  MATMLT,  TRMMLT , MULTTR, 

MATVCA , AND  LINSLV 

The  subroutines  MATMLT,  TRMMLT,  MULTTR,  MATVCA,  LINSLV,  and  the 
following  recursion  relationships  are  taken  from  a research  report  by 
D.  H.  Sinnott  [5] . These  subroutines  represent  an  efficient  algorithm 
for  solving  equation  (1-15) 

[ywg  + YhS]^  = +imp 

W2  tlS 

where  [Y  6 + Y ] is  a symmetric  block-Toeplitz  admittance  matrix, 
two  main  advantages  of  this  algorithm  are: 

(1)  It  solves  equation  (1-15)  directly  without  inversion 
of  [YWg  + Y^S]  and  is,  therefore,  more  efficient. 

(2)  Storage  requirements  are  considerably  less  than  that 

required  for  an  inversion  solution.  Approximately 

W2.  VlS  ”1 

one  quarter  of  the  matrix  [Y  + Y ] must  be  stored 
for  an  inversion  solution  while  only  the  first  row  of 
blocks  of  [YW®  + Y^S]  is  required  for  the  given  algo- 
rithm solution. 

Recall  the  form  of  the  block-Toeplitz  matrix  Y, 


[Y]  = Y(n)  = 


where  Y^  is  the  submatrix  which  defines  the  self-admittance  of  an  element 

of  the  array  (specified  number  of  apertures  - see  Fig.  1-5)  and  Y | ^ ^ j , 

i ^ j , is  the  submatrix  which  defines  the  mutual  admittance  between  ele- 

( s ) 

ments  i and  j of  the  array.  The  parenthesized  superscripts,  as  in  Y , 

are  used  to  identify  the  order  of  the  matrix  since  the  algorithm  to  be 

presented  is  defined  recursively  on  n.  Since  Y^  is  also  symmetric,  all 

submatrices  Y1  are  symmetric.  The  submatrices  are  N * N complex  matrices 
i P P 

where  N is  the  shortest  row  or  column  dimension  of  the  array  lattice. 


Y Y 

0 1 

Y Y 

1 0 


n-1 


Y Y ...  Y 
n n-1  0 


(1-108) 


(1-15) 

The 


Since  the  submatrices  are  of  order  N x N , then  Y is  of  order 

P P 

N • (n+1)  x N • (n+1).  A recursive  system  of  equations  defines 
P P 

further  sets  of  N x N matrices, 

P P 


^(0)  = Y_1  Y 
o o 1 


A(-X)  = Y-1 
o 


a<”-1)  - ti„  - (*(m:1))2]"1  4<"-2) 

N m-1  J 


(m)  _ a (m-1) 


m-1 


C'  = - A'"'  i l * Y 


(m-l)T 


s=0 


- Y ] 
m-s  m+1 


(m)  = (m-1)  _ (m-1)  (m) 

r rr  m-r-1  m 


0 < r < m-1 


J 


for  m=l,2, . . . ,n-l  (superscript  T denotes  transpose).  Then  the  N x n 


-1 


submatrices  of  Z = Y , defined  by 


are  given  by 


II 

o 

o 

N 

A(n-i) 

( 

II 

o 

u 

N 

- 

r-l 

A(n-1) ^ i < r < n 

( 

Z 

Z , 

+ A(n-1}  *(n:1)T 

rs 

r-l,s-l 

r-l  s-1 

. (n-1)  . (n-l)T  . „ 

A 1 < r,  s < n 

( 
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where  ^ = - I , the  unit  matrix  of  order  N and  tl/n  ^ = 0 
-i  N p r 

for  r ^ n.  ^ 

■f 

The  solution  for  V,  partitioned  into  subvectors  of  length  N 

F 

denoted  V (r=0, 1 , . . . ,n)  is 


vr  = I z 

L r\ 


s=0 


rs  s 


where  I is  the  sth  subvector  of  I of  length  N . Next  define  the 
s P 

matrices  <b  where 
r 

4>r  = A(n_1)  ^n~1)T,  r-0,1 n-1 

and  the  vectors  a and  b where 
r r 


n-r 


1 = l <f>  , I , 

r s-1  s+r 

s=0 


r r-0,1... 


n-r 

b = i 

r s=0 


*n-s-l  s+r 


,n 


so  that  the  general  solution  for  r > 1 is 


r r 

v = y ^ a - y ^ b , , . 

r rs-l  r-s  Tn-s  r-s+1 

s=0  s=l 

The  algorithm  can  be  briefly  summarized  by  the  following 


steps : 


(1)  Use  the  recursion  equations  (16)  to  define  tj/n  ^ 

.(n-1)  r 

and  A 

(2)  Form  the  matrices  <p^  by  equation  (22). 

(3)  Form  the  vectors  a^  and  b^  by  equation  (23). 

(A)  Form  the  voltage  subvectors,  V , by  equation  (24) . 


(21) 


(22) 


(23) 


(24) 
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Using  this  algorithm,  the  number  of  multiplications  and  divisions 

required  to  invert  Y^n^  varies  as  (n+l)“  • with  increasing  n and 

N.  This  should  be  compared  with  a variation  as  [(n+1)  • N for  a 

P 

general  elimination  method  and  demonstrates  a considerable  improve- 
ment when  n is  large. 

Since  the  computer  program  was  taken  from  [5],  no  further 
description  of  it  will  be  provided  other  than  mentioning  that  there 
are  numerous  comment  cards  in  LINSLV  written  using  the  same  notation 
as  the  algorithm  just  described.  For  the  subroutines  MATMLT,  TRMMLT, 
MULTTR,  and  MATVCA,  comment  cards  which  describe  the  matrix  multipli- 
cation operation  performed  follow  the  subroutine  statement  cards. 

In  summary,  the  subroutine  LINSLV(CE,  V,  YS,  NP,  NW)  solves  equa- 
tion (1-15)  for  V given  [YW^  + y^s]  and  I1™*3.  The  arguments  of  the 
subroutine  statement  card  represent  the  following  parameters:  the  first 

N elements  of  CE  are  elements  of  Iim^,  V is  the  output  magnetic  current 

we  hs 

coefficient  column  matrix,  YS  is  the  input  matrix  [Y  + Y ] partitioned 
in  terms  of  blocks  of  dimension  NP  x NP  (only  one  row  of  blocks  are 
required  [Y^,  , . . . , Y^]),  NP  is  the  smaller  of  the  row  or  column  dimen- 

sions for  the  array  lattice,  and  NW  uses  the  same  definition  as  NP  but 
replaces  the  word  smaller  with  larger. 

Minimum  allocations  are  given  by 

COMPLEX  A(NP*NP),  B(NP*NP),  C(NP*NP) 
in  subroutines  MATMLT,  TRMMLT,  MULTTR,  and  MATVCA,  and  by 
COMPLEX  CE(N+NP*NW+L) , PS ( 2*(NW+1 ) *NP*NP) , 

V (N) , YS ( (NW+1 ) *NP*NP) 
in  subroutine  LINSLV. 


C LISTINGS  OF  THE  SUBROUTINES  MATMLT • TRMMLT • MULTTR • 

C MATVCA.  AND  LINSLV 

C 

SUBROUTINE  MATMLT (A. B.C.NPI 
C CALCULATES  C=A*B 

COMPLEX  A<4).B<4)  .C(  4)  .D 

I J=0 

L=1 

DO  10  1=1.  NP 
DO  1 I J=1.NP 
I J=  I J«-  1 

o= ( 0 . . 0 . I 

KJ=L 

JK=J 

00  12  K= 1 * NP 
D=D+ AC JKJ*B( K J) 

JK=  J K'fNP 

KJ=K  J*1 
12  CONTINUE 
C(  I J )=D 
It  CONTINUE 
L=L+NP 

10  CONTINUE 
RETURN 
END 

SUBROUTINE  TR MMLT ( A . B. C . NP) 

C CALCULATES  C=- TRANSPOSE! A) * B 

COMPLEX  Al 4 ) ,B<4) ,C! 4) .D 

1 J = 0 
L=1 

DO  10  1=1. NP 
M=1 

DO  11  J=1.NP 
IJ=I  JM 
D= ( 0 . . 0.  ) 

K I = M 
KJ=L 

DO  12  K=1.NP 
D=D*A( KI )*b(KJ ) 

KI=KH-1 
KJ=K  J ♦ 1 
12  CONTINUE 
M=M*NP 

C! I J ) =C ( IJ ) -D 

11  CONTINUE 
L=L4NP 

10  CONT 1 NUE 
RETURN 
END 

SUBROUTINE  MULTTRI A, B. C. NP» 

CALCULATES  C = A*TR ANSPO SE  < B) 

COMPLEX  AI4).B<4) ,C<4» .0 
I J=0 

DO  10  1=1. NP 
DO  10  J= 1 . NP 
I J=I J4 1 
D=(0.,0.) 

IK=I 
JK=J 

00  11  K=1.NP 


C 


n n 


0*0 A(  JK ) *U  < IK) 

IK=IK*NP 
JK=JK4-NP 
11  CONTINUE 
CC1 J)=D 
10  CONTINUE 
RETURN 
END 

SUBROUTINE  MAT VCA I A , B. C . N » N 1 1 


POSITIVE  ACCUMULATION  OF  A *B  IN  C IF  Nl=l  OR  NEGATIVE  ACCUMULATION 
A*B  IN  C IF  Nl=2 
COMPLEX  A( 4 ) .814)  .C(  4)  .0 
00  10  1=1 . N 


0=10. .0.1 
I J=  I 

DO  11  J= 1 . N 
D=D*A(  I J )*B(  J) 

I J=  I J-f  N 

11  CONTINUE 

GO  TO  (12.13I.N1 

12  C( I >=C( I >«-D 
GO  TO  10 

13  CU  >=C(  D-D 
10  CONTINUE 

RETURN 

END 

SUBROUTINE  L I NSL V CCE • V • VS. NP ,NM) 

COMPLEX  CE ( 13>.PSI32).VI6). VSf 16) 

01  ME  NS  I ON  I A ( 2 ) 

100  FORMAT1/V1 OX. • ILLEGAL  CALL  TO  LI NSL V — — NK  = »»I4> 

101  FORMAT (//l OX. • ILLEGAL  CALL  TO  LINSLV  - - NP  =»,I4) 

IFINW.LT. 2)  GO  TO  10 

IF(NP.LT.2)  GO  TO  11 
C CALC  DEL  { — 1 * AND  PSUOI.O) 

N=NW-1 
NPW=NP*NW 
N2=NP*NP 
DO  12  1=1. N2 
PS1 1 )= YS 1 I ) 

12  CONTINUE 

CALL  LINEQ(PS.NP) 

CALL  MATMLT (PS. YS1N241 ).PS(N2+1). NP ) 

C IA( 1 ) = START  ADDRESS  IN  PS  ARRAY  OF  PSl(M-l).O) 

C I A( 2 ) = START  ADDRESS  IN  PS  ARRAY  OF  PS((M).0) 

C IA11I-N2  = START  ADDRESS  IN  PS  ARRAY  OF  DEL1M-2) 

C IA121-N 2 = START  ADDRESS  IN  PS  ARRAY  OF  DEL 1 M— 1 ) 

C ISTFi=2*NW*N2Fl  = START  ADDRESS  OF  AN  ADDITIONAL  SCRATCH  AREA 

I A( 1 )=N2*1 
I A ( 2 )=NW*N2 +N2F1 
I ST=2*NW*N2 
MZ=N2«-N2M 
MM=0 

C ITERATE  ON  M=1.2....N.  FOR  M*N.  ONLY  CALC  DEL  < M — l ) . 

DO  13  M=  1 . N 
10=1  A(  1 )*MM 
MM=MMFN2 
I 1= I A( 2 ) ♦MM 

C 10  IS  START  ADDRESS  OF  PSC I M-I ) . M- 1 ) 

C II  IS  START  ADDRESS  OF  PSl(M).M) 

C CALC  OEL(M-l) 


! 
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CALL  MATMLT (PSt 10). PSt  I 0 I . P SC  I ST ♦ 1 ) . NP I 
IJ=I ST 

DO  14  1=1. NP 

00  14  J=1.NP 
I J=I J*1 

PS« I J>=-PSt 1 J) 

IFtl.EQ.J)  PS(IJ)=PS(IJ)M. 

14  CONTINUE 

CALL  L I NEQ ( PSC 1ST ♦ 1 J .NP) 

10=1 A{ 1 »-N2 
ID1=IA«2)-N2 

CALL  MATMLT (PSt I ST*1 ) . PS ( ID  I *PS(  ID  1 1 . NP ) 

IFtM.EQ.N)  GO  TO  15 
C CALC  PS  t t M ) , M ) 

MZZ=  MZ 
MS=I A( 1 ) 

1 J=  I ST 

OO  16  I=I.N2 
I J=I J* 1 

PSt I J ) =YS( MZZ) 

MZZ=MZZ* 1 

16  CONTINUE 
MZZ=MZ— N2 

DO  17  I S= 1 * M 

CALL  TRMMLT (PSIMS).YSI  MZZI.PSI  I ST*  I) , NP ) 

MS=MS*N2 

MZZ=MZZ-N2 

17  CONTINUE 
MZ=MZ*N2 

CALL  MATMLT  t PSt  ID  1)  .PSt  I ST*  l l.PSIHI.NP) 

C CALCULATE  PSttMI.R)  FOR  R=0.l....  M-l.  (IR=R| 

I 0R= IAl  1 ) 

1 1 R=  I A ( 2 ) 

IMR=  10 

OO  18  I R= 1 . M 

CALL  MATMLT f PS t I MR ) . PS  t I 1 ) • PS t I ST *1 ) . NP I 
I MR= I MR— N2 
I J=I ST 

00  18  1=1. N2 

PSt  I 1R  )=PS(  I OR  )-PS(  I J*  1 ) 

IJ=I J*l 
1 1 R=  1 1 R* 1 

1 0R=  I OR* 1 

18  CONTINUE 

I = I A t 1 » 

I At  1 )=  I A 1 2 ) 

I At 2 )=I 
13  CONTINUE 

C HAVE  FINISHED  ITERATION  ON  PS.  NOW  PUT  PHItR)  INTO  PStIAt2)) 

15  I PH I =1 At2)-N2 
I PSI  = I At  II 
IOR=  lAt 1 I 

I IR=  1 At  21 
OO  19  1=1.  N 

CALL  MUL  TT  R t PS  t I PHI  I. PSt  I OR).  PSt  URl.NP) 

I 0R= I 0R*N2 
1 1R=  I 1 R*N2 

19  CONTINUE 

PUT  PH  I t — 1 ) IN  PS  t I PHI  ) 

J= IPHI 


C 


r\  o 


-J*'1  ■ ' WW  * 


c 


DO  20  I=1,N2 
PSCJ )=-PSI  J) 

J=J*  1 

20  CONTINUE 

NOW  HAVE  PHI IS),S=-1 .0.1 . .. • N— I STARTING  AT  PS(  IPHI  I 
AND  PSIIS).S=  0*1. 2*. ..  N-l  STARTING  AT  PSIIPSII 

IB=NW*NP+1 
10=1 

J=2*NW*NP 
DO  21  I* 1 • J 
YS( I >=(0.,0.) 

21  CONTINUE 
IV=1 

DO  22  J=1.NW 
NR— N W—  J* l 
IIS= IPHI 
1 2S=  IPHI  4-N4N2 
IVS=IV 

DO  23  1=1, NR 

CAUL  MATVCAIPSII1S),  CE ( IVS) > YS(  IC),NP.l  I 
CALL  MAT VCAI PSI I2S).CE( IVS) * YSI IB),NP.I) 

1 1 S= I 1 S*N2 
I2S=I2S-N2 
IVS= IVS ♦NP 

23  CONTINUE 
IB=IB«-NP 
IC=IC*NP 
I V=I V+NP 

22  CONT  INUE 

NOW  CALCULATE  V IN  I LOCATIONS 
J=NW*NP 
DO  24  1=1  . J 
CE ( I ) =—  Y S ( I ) 

24  CONTINUE 
I V=NP  + I 

I B= I V+J 
IC=1 

DO  25  I R=1 . N 
I 1S=IPSI 

I2S=IN-i )*N2*IPSI 
I BS= I B 
ICS=  IC 

DO  26  IS=1  . IR 

CALL  M ATVCA  < PSI  I 1 S) • YSI ICS)  • CE I IV) .NP.  I ) 

CALL  MAT VCAI PS I I2S) . YSI IBS) .CEI I V) .NP.2) 

I 1S=  II  S*N2 
1 2S=  I2S-N2 
ies= I8S-NP 
ICS=  ICS-NP 

26  CONT INUE 
I B=  I B*-NP 
IC=IC*NP 
IV=  I V4-NP 

25  CONTINUE 

DO  27  1=1, NPW 
VII )=CEt I ) 

27  CONTINUE 
RETURN 

10  WRITEI3.100)  NW 
RETURN 

11  WRITE(3,101)  NP 
RETURN 

END 
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IX.  DESCRIPTION  OF  THE  SUBROUTINES  LINEQ,  DECOMP,  AND  SOLVE 

The  subroutine  LINEQ(C,LL)  is  used  in  subroutine  LINSLV  to  invert 
a complex  matrix.  The  input  to  LINEQ  consists  of  a square  complex 
matrix  of  order  LL  x LL  stored  columnwise  in  C and  the  dimension  vari- 
able LL.  The  output  from  LINEQ  is  C ^ stored  columnwise  in  C. 

Minimum  allocations  are  given  by 

COMPLEX  C(LL*LL) 

DIMENSION  LR (LL) . 

The  subroutines  DECOMP (N , IPS ,UL)  and  SOLVE (N , IPS ,UL,B ,X)  are 
called  from  the  main  program  to  solve  the  matrix  equation 

[YW8  + Yhs]V  = tlmP  (1-15) 

when  the  array  lattice  is  isosceles  triangular . This  subroutine  combi- 
nation uses  the  method  of  Gaussian  elimination  and  LU  decomposition 
described  in  [6,  Section  9] . The  input  to  DECOMP  consists  of  N and  the 

\jp  hs 

N by  N matrix  of  coefficients  [Y  + Y ] which  is  stored  by  columns  in 

UL.  The  output  from  DECOMP  is  IPS  and  UL.  The  output  is  fed  into  SOLVE. 

The  rest  of  the  input  to  SOLVE  consists  of  N and  the  column  of  coef- 
-*imp 

ficients  I stored  in  B.  Solve  puts  the  solution  V to  the  matrix  equa- 
tion in  X. 

Minimum  allocations  are  given  by 

COMPLEX  UL(N*N) 

DIMENSION  SCL(N) , IPS(N) 

in  subroutine  DECOMP , and  by 

COMPLEX  UL(N*N) , B(N),  X(N) 

DIMENSION  IPS(N) 

in  subroutine  SOLVE. 

3 

DECOMP  and  SOLVE  require  roughly  N / 3 multiplications  and 

to  solve  a system  of  N linear  equations  whereas  using  LINEQ  in 

3 2 

which  requires  an  inverse  needs  roughly  N + N multiplications 
divisions  to  solve  the  same  system  of  equations. 


divisions 
a method 
and 
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SUBROUTINE  LINEQ(C.LL) 

COMPLEX  C(  4) .STOR.STO. ST.S 
DIMENSION  LR(2) 

DU  10  I=1.LL 
LR ( I » = I 

10  CONTINUE 
M 1 = 0 

DO  11  M= 1 * L L 
K=M 

K2  = M 14-K 

51  = ABS(REAL(C(K2)  ))4-ABS(  AIMAG(C(K2>)  ) 
DO  12  I=M,LL 

K1  = M1H 

52  = ABS  C REAL  (C(  K1  ) ) ) 4- AB  S ( A I M AG  ( C ( K 1 ) ) ) 
IFIS2-SD  12.12.13 

13  K=I 
S1  = S2 

12  CONTINUE 
LS=LR( M) 

LR(M)=LR(K) 

LR(K  ) = L S 
K2=  M 1 4-K 
S TOR  =C C K2 J 
J1  = 0 

OO  14  J= 1 • LL 
K1  = J 1*K 
K2=  J I 4-M 
STO=C ( K 1 ) 

C(K1  )=C(K2  ) 

C(K2)=ST0/ST0R 
Jl=J 14-LL 

14  CONTINUE 
K 1=M  14-M 

C < K 1 )=1  ./ST  OR 
DO  15  1=1. LL 
IF!  I — M ) 16.  15,16 

16  K 1 =M  1 4- 1 
ST=C I K 1 ) 

C I K 1 ) = 0 • 

J1  = 0 

DO  17  J= 1 • LL 
K 1 = J 1 ♦ I 
K2-J  14-M 

CIKl )=C ( K 1 >-C(K2»*ST 
Jl  = J 1 4-LL 

17  CONTINUE 

1 5 CONT INUE 
M1=M14-LL 

11  CONTINUE 
J1  = 0 

DO  18  J=1.LL 
IF ( J — L R ( J ) ) 19.20.19 

19  LRJ=LR ( J ) 

J2=(  LR J- 1 ) *LL 
DO  21  1=1. LL 

K2= J24- 1 
Ki  = J 14-1 
S=C( K2 ) 
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CCK2 J=C(K1  ) 

CtKl  ) = S 
21  CONTINUE 

LR<  J ) =LR(LR  J ) 

LR(LRJ)=LR J 
IF { J— LR ( J) ) 19.20*19 

20  Jl  = JI+Lt_ 

18  CONTINUE 
RETURN 
END 

SUBROUTINE  DECOMP ( N. IPS. UL) 

COMPLEX  UL ( 36) .PI VOT .EM 
DIMENSION  SCLI6) . IPS (6 » 

DO  5 I=1.N 
IPS< I ) = I 
RN=0  - 
J1  = I 

DO  2 J=1 »N 

ULM  = ABS< REAL! UL(  J 1 ) ) > + ABS < AI MAG < UL < J 1 ) ) 1 
J1=J 1 ♦N 

IF(RN-ULM)  1.2.2 

1 RN=ULM 

2 CONT  INUE 

SCL ( I >=1./RN 
5 CONT  INUE 
N M 1 = N—  1 
K2=0 

DO  17  K=l,NMl 
8IG=0. 

DO  11  I=K,N 
IP=I PS ( I ) 

I PK—  IP4-K2 

S I ZE  = ( ABSi  REAL  ( UL  ( IPK)  ) ) 4-AB  S ( A I MAG  ( UL  ( I PK  ) ) ) ) * SCL  ( I P ) 
IF(SIZE-BIG)  11*11.10 

10  BIG=SIZE 
IPV=  I 

11  CONTINUE 

IF 1 I PV—K ) 14.15.14 

14  J=IPS(K) 

IPS! K)  = I PS ( IPV) 

IPSC IPV)=J 

15  KPP— I PS  1 K ) ♦ K2 
PI VOT=UL{ KPP ) 

KP1=K+1 

DO  1 6 I=KP1  .N 
KP=KPP 

IP=IPS<  I ) + K2 
EM=-UL( IPJ/PIVOT 
18  UL  (IP)  =— EM 

DO  1 6 J=KP1 .N 
IPs I P+N 
KP=KP+N 

UL ( IP! =UL( IPI+EM4UL1KP1 

16  CONTINUE 
K2=K24N 

17  CONTINUE 
RETURN 
END 

SUBROUTINE  SOL VE< N . I PS .UL . B • X » 

COMPLEX  UL ( 36  1.6(61. X(  6). SUM 


DIMENSION  IPSI6) 

NPl=N+l 
IP=I PSI l I 
XI  1)=BI IP) 

DO  2 1=2. N 

IP=IPS( I ) 

IPB= IP 
I Ml  = 1 — 1 
SUM=0. 

DO  1 J= 1 . IM  1 
SUM=SUM+UL< IP)*X( J) 

1 IP=  I P*-N 

2 XI I)=B( IPB)-SUM 
K2=N* ( N— 1 ) 

IP= I PS ( N ) ♦K  2 
X(N)=XCN)/UL< IP) 

DO  4 IBACK=2.N 

I =NP 1 -I  BACK 

K 2=  K 2—  N 

IPI=  IPS! I ) +K2 

IP1=I+1 

SUM=0. 

IP=IPI 

DO  3 J= I P X . N 
IP=I  Pt-N 

3 SUM=SUM*UL( IP)*X( J) 

4 XI  I )=(  XI  I )-SUM)/UL<  IPI  ) 
RETURN 

END 
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X.  DESCRIPTION  OF  THE  SUBROUTINE  TOPGEN 

The  subroutine  TOPGEN (N,  NP,  NT,  NW,  X,  Y,  T,  Tl,  KR)  generates 
a complete  matrix  T given  one  column  Tl. 

Minimum  allocations  are  given  by 

COMPLEX  T(N*N) , T1(N),  T2(N,N) 

DIMENSION  D(6) , X(6) , Y(6) 

The  variables  in  the  argument  of  TOPGEN  are  defined  as  follows: 

N is  the  dimension  of  the  square  matrix  T,  N is  the  shortest  row  or 

P 

column  dimension  of  the  array  lattice,  NT  specifies  the  array  lattice 

type-1  for  rectangular  and  2 for  isosceles  triangular,  NW  is  the 

largest  row  or  column  dimension  of  the  array  lattice,  X(I)  and  Y(I)  are 

the  center  coordinates  of  the  Ith  aperture  and  T is  the  output  matrix 

which  is  generated  from  one  column  Tl  of  T.  KR  is  a parameter  which 

2 2 

determines  whether  or  not  D(I)  is  calculated  (D(I)  = (X^-X^.)  + (Y  -Y^)  ) . 
If  KR  = 1,  D( I)  is  evaluated  whereas  if  KR  = 2,  this  step  is  skipped 
(D(I)  has  been  determined  already  by  a previous  call  to  TOPGEN) . 

DO  loops  11,  12,  13,  and  14  generate  [T]  given  tT  for  a uniformly 
spaced  rectangular  array  lattice  while  DO  loops  16,  17,  18,  19,  21,  and  22 
does  the  same  for  a uniformly  spaced  isosceles  triangular  array  lattice. 
For  the  rectangular  array  lattice,  [T]  which  is  symmetric  block-Toeplitz 
is  determined  by  forming  one  column  of  blocks  and  generating  the  rest  of 
[T]  using  the  Toeplitz  property.  For  the  isosceles  triangular  array 
lattice,  the  elements  T„  below  the  diagonal  of  [T]  (excluding  the  ele- 
ments of  the  first  column  which  are  given  (Tl))  are  determined  by  com- 
paring the  distance  squared  between  the  centers  of  the  apertures  i and  j, 

2 2 

D1  (D1  ofT..  “ (x.-x.)  + (y.-y.)  ),  to  the  distance  values  associated 

J r j 1 * - 

1_  , . C *.1.  _ f J . . _ 1 Tv  / Y-v  /“  m / \ ^ . / \^\ 


with  the  elements  of  the  first  column,  D(D  of  T..  = (x,-x.)  + (y, -y.)  ). 

il  1 i W1 

When  D1  = D,  the  value  used  for  T^  which  has  an  associated  D1  value  is 
the  same  as  the  corresponding  element  of  the  first  column  with  an  associ- 
ated D value.  The  elements  above  the  diagonal  are  determined  from 
symmetry,  y = y (consequence  of  Galerkin's  method). 


SUBROUTINE  TOPGEN IN. NP . NT • NW • X • Y .T . T I . KR ) 

COMPLEX  T< 361 .Til  6) .T216.6) 

DIMENSION  D(6).X(6).Y< 6) 

npw=np*nw 

NPKPcNPW*NP 
N1=N— 1 

IF1NT.EQ.2)  GO  TO  10 
DO  1 1 1=1 • NW 

L 1 = 1 

DO  12  J=  1 ♦ NP 
DO  12  K=1.NP 
L = 1 + IABS1K— J> +NP* 1 1 — 1 ) 

T2<  I ,L1  >=T1  (L) 

Ll=Ll+l 

12  CONTINUE 
1 1 CONT INUE 

DO  13  1=1. NW 
DO  13  I 1=1 .NW 
1FII.EQ.1)  11=11 
IF1I.GT.1)  I 1=1-1 1+1 

IF(  I l.GT.l  ) .AND.  I II.LE  .0)  ) 11=11-14-1 

12=1 

DO  14  J= 1 . NP 
DO  14  K=1.NP 

L 1 = 1 I-l)*NPWP+<  U-l)*NP+l  J-l  ) 4NPW4K 
T(L1 »=T2(1 1.12) 

12=1 2+1 

14  CONTINUE 

13  CONTINUE 
RETURN 

10  IF ( KR  .EQ  .2 ) GO  TO  IS 
DO  16  1=1, N 

D<  I )=(  XI  I) -XI  1 J )*«X<  I)— XI  1 ) )♦<  VI  D— r<  1 ) >4<Yl  I )-Y<  1)1 

16  CONT  INUE 

15  11=1 
K 1=  1 

DO  17  1=1. N1 
DO  18  J= I 1 * N 

01=1X1 J)-X<  XI 1*1 XI Jl-Xl  I))  + (Y1 J>— Y1 I))*1Y< J)-Y(  I)> 

DO  19  K=1.N 

IF1ABS1D1— D(KI).LT.O.Ol)  T(K1)=TKK) 

IF ( ABSl Dl-Dl K ) ) .LT.O .0 1 ) GO  TO  20 

19  CONTINUE 

20  K1=K l+l 
18  CONT INUE 

IF(I.EQ.l)  11=11+1 
11=11+1 
K1=K1+I 1-1 

17  CONTINUE 
J2=0 

DO  21  J= 1 , N 
J 1 = J 

DO  22  1=1 • J 
J3*J2+I 

IF ( I-J ) 23. 24.23 
24  T!  J 1 )=T1  ( 1 ) 

GO  TO  22 
23  T { J3 )=T 1 J1 I 
J 1= J 1 + N 
22  CONT INUE 
J2= J2+N 

21  CONTINUE 
RETURN 
END 
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XI.  DESCRIPTION  OF  THE  MAIN  PROGRAM  WITH  SAMPLE  INPUT-OUTPUT  DATA 

The  main  program  computes  the  x,y  coordinates  of  the  centers  of 
the  apertures  in  the  array,  the  complex  coefficients  which  deter- 
mine the  magnetic  currents  M1  according  to  (1-6)  and  the  scattering 
coefficients  S according  to  (1-103)  and  (1-104).  The  main  program 
calls  the  subroutines  AY,  YHSP,  TOPGEN,  CSMTZ,  DECOMP,  SOLVE,  and  LINSLV. 

The  data  cards  are  read  into  the  main  program  according  to 

READ  (1,100)  N,  AL,  AL1 , BL,  BL1 , ER,  LM,  LN 
100  FORMAT  (14,  5F7.4,  213) 

READ  (1,102)  NC,  NR,  NT,  NE,  DX,  DY , DT 
102  FORMAT  (414,  3F7.4). 

The  variables  AL,  AL1,  BL,  and  BL1  are  respectively  a/X,  a'/X,  b/X,  and 

b'/X  where  X is  the  free  space  wavelength.  ER  is  the  relative  dielectric 

constant  c of  (1-34)  and  (1-35)  inside  the  waveguide.  The  variables  LM 

and  LN  are  respectively  the  total  number  of  m and  n modes  used  in  deter- 
TE  TM 

mining  A^  (1-30)  and  A^  (1-31)  (only  odd  m starting  with  m=l  and  even 

n starting  with  n=0  (TE)  or  n=2  (TM)  are  considered  due  to  the 

sin  ^ cos  factor  appearing  in  (1-33)).  LM  which  represents  the  con- 

W2 

tribution  of  the  mth  waveguide  mode  to  Y should  be  chosen  so  that  the 

2 2 2 

contribution  of  the  (l/(m  /a  - 1/a'  ))  cos  (mTTa'/2a)  factor  in  (1-33) 

TE  TM 

results  in  very  small  A^  and  A.^  values.  LN  which  represents  the  con- 
tribution of  the  nth  mode  to  YW^  should  be  chosen  so  that  the  argument 
n7rb'/2b  of  the  sin(  )/(  ) factor  in  (1-33)  is  greater  than  7r.  NC  is  the 
number  of  columns  of  apertures  measured  in  the  x direction.  NR  is  the 
number  of  rows  of  apertures  measured  in  the  y direction.  NT  specifies 
the  array  lattice  type,  1 for  rectangular  and  2 for  isosceles  tri- 
angular. NE  is  the  driven  aperture  number.  DX  is  the  distance  per  unit 
wavelength  between  the  closest  outer  waveguide  edges  in  the  x direction 
while  DY  Is  the  same  but  in  the  y direction  (see  Fig.  1-1).  DT  is  the 
waveguide  wall  thickness  per  unit  wavelength. 
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Minimum  allocations  are  given  by 

COMPLEX  CE(N+NP*NW+1) , S(N*N),  S3(N), 

V(N),  Y0(2*LM*LN+LM) , YS(N*N), 

YS1 ( (NW+1) *NP*NP) , YSP(N) 

DIMENSION  A( 2 *LM*LN) , IPS(N),  X(N) , Y (N) 

in  the  main  program,  by 

COMPLEX  Y(2*LM*LN  + LM) 

DIMENSION  A(2*LM*LN) , F(LM) , STE(LN) , STM(LN) 
in  the  subroutine  AY,  by 
COMPLEX  YSP(N) 

DIMENSION  X(N),  Y(N) 
in  the  subroutine  YHSP,  by 

COMPLEX  T(N*N),  T1(N),  T2(N,N) 

DIMENSION  D(N) , X(N) , Y (N) 
in  the  subroutine  TOPGEN,  by 

COMPLEX  A(N),  B(N),  E(N) , ES(N-l),  X(N) 
in  the  subroutine  CSMTZ,  by 
COMPLEX  UL(N*N) 

DIMENSION  SCL(N),  IPS(N) 
in  the  subroutine  DECOMP,  by 

COMPLEX  UL(N*N) , B(N) , X(N) 

DIMENSION  IPS(N) 
in  the  subroutine  SOLVE,  and  by 

COMPLEX  CE(N+NP*NW+1) , PS(2*(NW+1)*NP*NP) , 

V (N) , YS((NW+1)*NP*NP) 

in  the  subroutine  LINSLV  where  NP  is  the  smaller  of  the  number  of  rows  or 
columns  of  apertures  for  the  rectangular  array  lattice  while  NW  uses 
the  same  definition  but  replaces  the  word  smaller  by  larger. 

DO  loop  11  calculates  the  x and  y coordinates  of  the  centers  of 
the  apertures  if  the  array  has  a rectangular  lattice  while  DO  loop  15 
does  the  same  if  the  lattice  is  isosceles  triangular.  Referring  to 
equations  (1-30)  to  (1-35),  statement  40  stores 
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A^  in  A(k  = (m+l)/2  + n/2  * LM) 

TF 

Y^  in  YO(k  = (m+l)/2  + n/2  * LM) 

m = 1,3,5,...,  LM 
n = 0,2,4, ... , LN 
TM 

Ak  in  A(k  = LM*LN  + (m+l)/2  + (n-2)/2  * LM) 
TM 

Yr  = Y0(k  = LM*LN  + (m+l)/2  + (n-2)/2  * LM) 

m = 1,3,5,...,  LM 
n = 2,4,6,...,  LN  . 


Note  that  only  odd  m and  even  n modes  are  calculated  and  that  the  first 
TE  TM 

subscript  of  A^  and  A^  in  (l-30)  and  (1-31)  has  been  dropped  since 

the  same  expansion  function  is  used  in  every  waveguide  region  - 

A^  = A2k  = = = DO  loop  18  calculates  YWg.  Statement  41 

hs 

stores  the  first  column  of  Y in  YSP.  Statement  42  generates  the  com- 
tis 

plete  Y matrix  given  the  first  column  and  stores  it  in  YS  by  columns. 
DO  loop  19  and  statement  43  generate  ilmp  defined  by  (1-106)  and  (1-107) 
and  store  it  in  CE. 


If  the  array  is  linear,  DO  loop  21  is  used  to  calculate  [Yw^  + Y^S] 
and  statement  44  calculates  the  coefficient  vector  V defined  by  (1-15) 
and  stores  it  in  V. 


If  the  array  contains  at  least  two  rows  and  columns  and  has  a 

rectangular  lattice,  statement  20  and  DO  loops  24,  25  calculate  the 

symmetric  submatrix  blocks  of  [YWg  + YhS]  given  [YWR]  and  a column  of 
Vis 

[Y  ] . If  the  array  contains  at  least  two  rows  and  columns  and  has  an 
isosceles  triangular  lattice,  DO  loop  27  calculates  [YWS  + Y*18]  given 
[YW^’]  and  [Y^S]  . Once  [YW^  + Y^S]  has  been  determined  for  these  two 
cases,  statement  26  stores  V in  V for  the  rectangular  array  while  state- 
ments 45  and  46  store  V in  V for  the  isosceles  triangular  array. 


n n 


42 


DO  loop  31  calculates  one  column  of  the  scattering  matrix  defined 
by  (1-103)  and  (1-104)  and  stores  it  in  S3.  Statement  47  generates  the 
complete  scattering  matrix  [S]  given  the  first  column  and  stores  it  in  S 

The  following  is  a listing  of  the  main  program  with  sample  input- 
output  data.  Figure  2 shows  the  coupled  power  (20  log  |S.^|)  and  phase 
of  S (i  = 2,6)  between  a driven  element  and  the  rest  of  the  elements 
for  the  2x3  waveguide-fed  rectangular  aperture  array  example  used  in 
the  sample  input-output  data. 


LISTING  OF  THE  MAIN  PROGRAM  AND  SAMPLE  DATA 

//XXXX  WAT  F I V (XXXX.XX.1.2)  ,'XXXX* .REGION=250K 
SJOB  XXXX • T I ME= 1 • PAGES— 30 

C 

C MAIN  PROGRAM 

C THIS  PROGRAM  CALLS  THE  SUBROUTINES  AY. YHSP. TOPGEN . CSMTZ . 

C DECOMP. SOLVE  * AND  LtNSLV 

COMPLEX  CE< 13) .G3.G5.G7.G9. S<36) .S3! 6) . U. VI 6) 

COMPLEX  YOI55)«YSt  36 I.YSIt 16).YSP<6).YWG 
DIMENSION  A(S0).IPS!6) *X(6).Y(6) 

CCMMCN  AL1 . 6LI .PI .PI  2. P 13, U/R I /ETA 

COMMON  /R3/G1 .G2.G3.G4 ,G5 «G6 ,G7 • G8  «G9 • G 10 • G1 I 

COMMON  /RS/AL .BL/R6/ER .LM.LN.NTS 

100  FORMAT! 14, 5F7. 4. 213) 

101  FORMAT !/5X, • N* • 5X . • AL« . SX . • AL I • . 4X . • BL • • 5X • • BL 1 • »4X, «ER» ,5X, *LM» • 
13X.*LN»/2X. I4.2X. 5F7.4 .2X, I3.2X. 13) 

102  F ORM  AT (414. 3F  7.4) 

103  FORMAT ( //4X . • NC * . 5X, • NR*  , 5X • • NT • • 5X, • NE • .SX.  «DX* . SX. • DY* *5X. »DT» / 
12X.I  4.3X , 14  *3X.  I4.3X.I 4.2X.3F7.4  ) 

104  FORM  AT ( /4X « •X*/(2X. 5EI 4.7)) 

105  FORMAT !/4X. •Y*/(2X»5E14.7)) 

106  F ORM  AT (/4X, • Y— WG* /2X.2E14.7) 

107  FORMAT !/4X, • Y- HALF  SPACE • /( 2X . 5E 1 4. 7) ) 

108  FORMAT !/4X. • V - - UNKNOWN  VECTOR • / ( 2X . 5EI 4 . 7 ) ) 

109  FORMAT! /4X. • S - - SCATTERING  MATRI X* / ( 2 X. 5E1 4.7)) 

READ ( 1 .100)  N. AL • AL I • BL. BL 1 * ER.LM.LN 

WRITE! 3.  10  1 ) N.AL.AL1 .BL.BL1 .ER.LM.LN 
READ  <1.1 02 ) NC.NR .NT .NE.DX.DY.DT 
WRITE! 3. 103  ) NC.NR. NT .NE.DX.DY.DT 
PI=3. 141593 
PI2=2.*PI 
P 1 3=  PI  /2  . 

ETA=376.730 
U=( 0 . . 1 .) 

N 1 — N ♦ 1 
N2=L  M4LN4 1 
NN=N*N 

AL2= AL/2.+DT 
BL2=BL/2.4DT 
DX2=OX/2. 


NP=A*«lN0<NC.NRI 
NW=A  MAX  0 ( NC  » NR ) 

NTS=2*L**LN 

NP2=NP*NP 

G1=AL1/PI 

G2=AL l /PI 2 

G3=-U*G2 

G4=G  1*G  1 

G5=-U*G* 

G6=G1*G2 
G 7-— U*G6 
G8=Gi  *G4 
G9=—  U*G8 

G 1 0= ( 1 • — 0.  25/ < AL l*AL I ) » 

G 1 1 = < 1 .♦0.2  5/C AL14ALI ) >*G1 

S1=DX*2.*DT+AL 

S2=DY *2 ■ *DT  +8L 

IF  C NT.EQ.2 ) GO  TO  10 

K=1 

OO  II  1=1. NW 
DO  1 1 J=1.NP 
IFCNC.GT.NR)  GO  TO  12 
X(K)=AL2+1 J-  1»*S1 
Y<K)=BL2*( I — 1 1 *S2 
GO  TO  13 

12  X(K)=AL2+C 1-1 )*S1 
Y ( K ) = SL2  ♦( J— l I *S2 

13  K=K«- 1 

11  CONTINUE 
GO  TO  14 
10  K=1 

DO  15  1=1. NR 
Kl=2-( 1-2*1 1/2 J 1 
DO  15  J=1.NC 

X(K)=C 14K1/2)*AL2+C J-l )*S1+K 1/2*DX2 
Y1K) =BL2+<  1-1  ) *S2 
K=K  + 1 

15  CONT  INUE 

14  WRITE13.104)  (X(I).I=1,N) 

WRITEC3. 105)  (YCI>.I=1.N) 

40  CALL  AYCA.YO) 

Y«G=f 0. ,0. ) 

K-  1 

DO  1 8 1=1 .NTS 
IFCI.EQ.N2)  K=K*LM 
YWG=YMG*A(  I ) *Y 0 ( K ) *A  C I I 
K=K>  1 

18  CONTINUE 
WRITEC3.106)  YMG 

41  CALL  YHSP(N.X.Y.YSP) 

KR=  1 

42  CALL  T OPGEN ( N . NP. NT • NW .X.Y.YS.YSP.KR) 
MR ITEC3.107)  IYSI I) .1=1. NN) 

00  19  1=1.  N 
CE ( I )=(0..0.) 

19  CONT INUE 

43  CE(NE)=2.*Y0C1 )*A(1 ) 

IFC < NC.GE.2) . ANO.CNR.GE.2) ) GO  TO  20 
YSPC  1 )*YSP<  1 > ♦ YWG 
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DO  21  1=1. N 

C£( 1 1 = CEI I 1/YSPI 1 ) 

IFC I .EQ.l ) GO  TO  21 
YSP  ( I ) - Y SP  ( I ly-YSPt  1 ) 

21  CONTINUE 

44  CALL  CSMTZI N.YSP.CE.V) 

GO  TO  22 

20  YSP< 1 )=YSP< 1 ) *Y*G 

IFINT.EQ.2)  GO  TO  23 

OO  24  1=1. NW 

DO  25  J=1.NP 

DO  25  K= 1 . NP 

L=( 1-1  I *NP2 ♦ ( J— 1 > *NP4K 

L 1 = 1 I- 1 >*NP+IABSCK- J1F 1 

YS1 ( L ) = YSP ( L 1 I 

25  CONT INUE 
24  CONTINUE 

26  CALL  LlNSLV(CE.V.YSl.NP.NW) 

GO  TO  22 

23  J=-N 

OO  27  1=1 ,N 
J=  JFN1 

YS( J>=YS< Jl+YWG 

27  CONTINUE 

45  CALL  DECOMPIN. IPS.YS) 

46  CALL  SOLVE ( N. IPS. YS.CE.V 1 

22  MR I TE( 3. 1061  (V(I).I=1,N) 

OO  31  1=1 . N 

S3! 1 ) = VI  I >*A<  1 ) 

IFII.EO.NE)  S3! I )=S3< I )— 1 . 

31  CONTINUE 
KR  *2 

47  CALL  T OP GEN I N < NP . NT . NW .X.Y.S.S3. KR ) 
MRITEI3.109)  (S(K),K=1.NN) 

STOP 

END 

SDATA 

6 1.0000  0.6500  0.4761  0.3095  1.0000  5 5 

3 2 1 1 0.9000  0.2000  0.0515 

tSTOP 
// 


PRINTED  OUTPUT 


N 

AL 

AL  1 

BL 

BL  1 

ER  LM 

LN 

6 

1.0000 

0.6500 

0. 4761 

0.3095 

1.0000  5 

5 

NC 

NR 

NT 

N£ 

DX 

DY  DT 

3 

2 

1 

1 

0.9000 

0.2000  0.0515 

X 

0. 551  5000E400  0 . 55 1 500 OE *■  00  0 • 2554499EF0  1 0 .2554499EF0  1 

0 . 4557498E+0 1 


V 

0 .269  5499E+0O  0 . 1 068649E+01  0 . 28 95499E4 00  0.  1 068649E f 0 1 

0. 1068649E+01 


Y-WG 

0. 1289030E-02-0.8829894E-04 


.4557498E+01 


.28954  99E+00 


J 


Y-HALF  SPACE 

0.  13751 10E- 02  0.6628551E-03-0. 34661 08E-03  0 . 1 4460 04E-0 3-0 . 1527445E-04 
-0.211 7092 E-0 5 0 . 6 1 5950 8E- 06  O. 1 8400 72E- 04- 0 . 374729 IE- 05-0. 1 479488E-06 
— O • 252 531 2 E— 05  0 .3 1 05473E- 05-0 .3466 1 08E-03  0 . 1 446004E- 0 3 0 . 13 751  1 OE- 02 
0.662 855 IE- 03  0 . 8 1 5950 8E- 06  0. 1 8400 72E- 04-0 . 1 52 7445E-04-0 .2117092E-05 
-0. 2525312E-05  0 . 3 1 05473E-05-0 .374 729 1 E- 0 5- 0 . 1 4 79488E- 06- 0. 1 527445E-04 
—0.211 7092 E— 05  0 . 8 1 5950 8E- 06  0 • 1 8400 72E- 04  0 . 1 375 1 1 OE- 02  0 .662855 IE- 03 
-0.3466 108E-03  0 . 1 44600 4E- 03-0 • 152 7445E- 04-0 . 21 1 7092E-0 5 0 . 81 595 08E-06 
0. I 84 0072 E- 04  0 . 8 1 5950 8E-06  0 . 1 840072E-04-0 . 1 527445E- 04-0 • 2 1 1 70 92E- 05 

— 0 • 3466 1 08E—  03  0 . 1 4460O4E- 03  0 • 1 375 1 1 0E-02  0 . 662 855 1 E-03  0 . 8 1 59508E-06 
0. 18400 72 E-04-0 . 1 52744 5E- 04-0 . 2 1 1 7092E-0 5- 0. 374 729 1 E- 05-0 . 1 4 79488E-06 

— O • 252 5312 E—  0 5 0 .3 1 05473E-05-0 . 1527445E-04-0 . 2 11 7092E- 05  0 .8 1 59508E-06 
0 . 1840072E-04  0 . 13751 1 OE- 02  0 .662 8551 E- 03-0. 3466 1 08E- 03  0 • 1 446004E-03 

— 0 • 252 53 1 2 E—  O 5 0 . 3 1 05473E-05-0 . 374729 1 E-05-0 . 14 79488E- 06  0 . 8 1 595 0 8E-06 
0. 1 84  00 72 E—  04— 0 . 1 52 744 5E- 04-0 . 2 1 1 7092E-05-0 . 3466 1 08E-0 3 0 . 1 4460 04E-03 
0.1375110 E—  02  0.6628551E  — 03 

I V - - UNKNOWN  VECTOR 

-0. 1237808E40 1 0 . 2902574E ♦ 00 -0 . 1171 50 IE* 00  0. 1 303046E+ 00-0 . 4 72 0550E- 02 
0 . 349 762 OE—  02  0 . 3 1646 1 2E-02  0 . 91 65 1 86E-02-0 . 1 3684 88E- 02  0 . 1 222063E- 02 
— 0. 577 7453 E— 03  0 . 2339229E- 02 

S - - SCATTERING  MATRIX 

-0. 7309630E— 01— 0 .21 7352 4E *00  0 .87725 1 0E-01 - 0. 97S7555E-01  0 .35348 73E- 02 
-0.2619111 E- 02-0. 2369 74 5E-02-0. 68631 32E-02  0 • 1 024760E- 02- 0 .9 1 5 1 1 3 1 E- 0 3 
0 . 432 631 OE- 03-0. 1751676E-02  0 . 87 7251 OE- 01 - 0. 9757555E- 0 1 -0 . 7309630E-0 1 
—0.21  73524  E*- 00— O .236974 5E— 02  — 0 .68631  32E-02  0 .35348 73E- 02-0 . 26  1 91  1 IE-02 
0.432631 OE-03-0.  1751676E-02  0 . I 024760E-02-0 . 9 1 5 1 1 31 E-03  0 .35348 73E-02 
-0.26191 1 IE-02-0 . 2369745E- 02-0. 6863 132E- 02-0. 7309630E- 0 1 -0 . 21 73524E400 
0.8772S1 OE—  01  — 0.9  75 7 55 SE— 01  0 .3534873E-02-0 . 26 1 91  1 IE-02-0 .2369745E-02 

— 0 . 6863 1 32E- 02— 0 . 2369745E  — 02  — 0 • 6863 1 32 E—  02  0 . 3534873E-02-0 . 26 1 9 1 1 IE-02 
0. 877251 OE— 0 1-0 . 9 75755 SE-Ol -0 . 7309630E-01-0 . 2 1 73524E ♦ 0 0-0 .2369745E- 02 

-O. 68631 32E- 02  0 . 35348 73E- 02-0.26 1 9 1 1 IE-02  0 . 1 024760E-02-0 .9 1 5 11 31 E-03 
0.432631 0E-03-0. 1751676E-02  0 • 3534873E- 02- 0. 26 1 91 1 IE—  02—0.23697  45E— 02 

— 0 • 6863 1 32  E— 02— 0 . 7309630E— O l — 0 .21 73524E*00  0 . 87725 1 OE- 0 1 - 0 . 97  57555E- 0 1 
0 .432 631 0E-03-0. 1 751676E-02  0 . 1 02 4 760E- 02-0 . 91 5 1 1 31 E- 03-0 .23697 45E-02 

-0 .6863132E-02  0 . 35348 73E- 02-0 -26 1 911 1 E-02  0 .8772510E- 0 1-0. 9757555E- 01 

I —0.7  309630E—  01  — 0.217  3524E  + 00 

; 


in  dB  while  the  lower  number  represents  the  phase  of 
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